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Abstract 


Fast Fourier Transform (FFT) and Least Square Error (LSE) estimation 
techniques were applied to the problem of identifying pilot-vehicle dynamic 
characteristics in flight simulation. A brief investigation of the effects of noise, 
input bandwidth and system delay upon the FFT and LSE techniques was 
undertaken using synthetic data. Data from a piloted simulation conducted at 
NASA Ames Research Center was then analyzed. The simulation was 
performed in the NASA Ames Research Center Variable Stability CH-47B 
helicopter operating in fixed-basis simulator mode. The piloting task consisted 
of maintaining the simulated vehicle over a moving hover pad whose motion 
was described by a random-appearing sum of sinusoids. The two test subjects 
used a head-down, color cathide ray tube (CRT) display for guidance and 
control information. Test configurations differed in the number of axes being 
controlled by the pilot (longitudinal only versus longitudinal and lateral), and in 
the presence or absence of an important display indicator called an 
acceleration ball". A number of different pilot-vehicle transfer functions were 
measured, and where appropriate, qualitatively compared with theoretical pilot- 
vehicle models. Some indirect evidence suggesting pursuit behaivor on the 
part of the test subjects is discussed. 
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1.0 Introduction 

A pilot when combined with a modern aircraft, whether an 
airplane or helicopter, forms one of the most complicated systems 
to be analyzed by a control engineer. Once the basic problem of 
stability is solved, the issue of handling qualities can be raised. In 
order to improve handling qualities, a better integration of man and 
machine is required, and in order to achieve better man machine 
integration, the system dynamics must be evaluated and optimized. 
The first step in this process is the measurement of the vehicle and 
pilot dynamics. The measurement of human dynamics is complicated 
by the fact that human characteristics are task dependent. Their 
dynamics can change dramatically depending upon the type of task 
being performed, and their familiarity with it. Dynamics can also 
vary between pilots for the same task. 

Early work in the area of human pilot dynamics (e.g. ref. 1) has 
led to the formation of a vast data base, which has aided many 
researchers in their development of models of the human pilot. Over 
the past three decades, many models have been proposed and tested. 
Models varying in complexity from the Crossover Model for single 
loop systems to the Structural Isomorphic Model. Although the 
crossover model is the simplest, it is the most general. When 
expressed mathematically the crossover model appears as, 


We*'* 

V Y - E — 
P T C ~ S 


“c 0 ' 18 

s 


( 1 ) 



2 


where Y p is the pilot transfer function, Y c is the plant transfer 
function, Kp and K c are gains, x is the system time delay, and co c is 

the crossover frequency. This model states that no matter what the 
plant dynamics, the operator compensates for them and the system 
crossover model is preserved. As the name implies this model is 
accurate near the crossover frequency, but if the low or high 
freqencies are important a different model must be used. The 
precision model for single loop systems (ref. 2) is accurate for a 
broader frequency range than the crossover model. If a more 
complicated model is required there are several multiloop examples 
to choose from. The McRuer Structural Isomorphic Model (ref. 2), the 
Linear Optimal Control Model (ref. 3), and the Hess Structural Model 
(ref. 4) are three options available. 

This paper uses simple single-loop models for generation of 
data that is used to test the least squares identification process and 
the Fast Fourier Transform analysis. Once the simple models were 
identified properly, a more complicated multiloop example was 
exercised. After the multiloop example was completed the 
identification procedure was used to identify vehicle pilot 
dynamics from data obtained from a CH-47B helicopter at the NASA 
Ames Research Center. 

2.0 Background 

2.1 Brief review of techniques for human transfer 
function estimation. 

When deciding upon an identification technique there are many 
choices, and depending on the system being identified there will be 
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advantages to using one method over another. Some of the different 
methods used employ orthogonal filters, spectral analysis 
techniques including Fast Fourier Transforms, and least squares 
estimation. 

The orthogonal filter method is a generalized technique that 
models the system dynamics as a series of transfer functions or 
linearly independent filters. The set of linerarly independent filters 
are of the form 


G(joo) - e 


-xjco h. 

j CO + 1 


P 2 (*,j<p- i) 

(TjO) + 1 ) ( X> jo) + 1 ) 




(2) 


where p r p 2 , . . . are determined by a regression technique, and T r 
. ■ ■ are predetermined time constants. The time constants are 

selected from models of the pilot that include the sensory organs, 
muscular mechanics, and the feedback therein, (see ref. 5 pg.4) 
Although general, the results from the orthogonal filter method are 
some what difficult to interpret due to the many parameters in the 
model. 

Power and cross power spectral densities can be computed to 
determine the pilot dynamics using spectral measurement 
techniques. The power spectral density of a random singal x(t) is 
derived from the autocorrelation function 

T 

«kxW 88 >oo 2T (3) 



which can also be defined as one half of a Fourier Transform pair 
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( 4 ) 


where <J> xx (cd) is refered to as the power spectral density of x(t), and 


®xx<“> = J <I> xx(' t ) e ' j “ t > 


( 5 ) 


Similarly the cross power spectral density of two random signals 
x(t) and y(t) is defined as 

a* . - 

J <I) xy(^ e "'* cot ^ : (6) 

•M 

Now if x(t) and y(t) are the input and output of a linear system, 
respectively, the transfer function of that system can be obtained 
as 


H(jco) ■ 


^ X y(Q>) 

<t xx(“) 



( 7 ) 


where H(s) is the system transfer function in the Laplace domain. 

If the input data can be described as a sum of sinusoids Eq. 7 
can be simplified and the transfer function can be obtained using a 
Fourier Coefficient method where 
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HGco) = 


Cy(jtQ) 

C x (jtD) 


( 8 ) 


where C y and C x represent Fourier coefficients whose real and 
imaginary parts are defined as 

i T i 

Re[C y (jcoj)] - y" Jy(t)sin(cojt)dt ( 9 ) 

'o 

i T i 

lnfiag[C y (jcoj)] - - Y“J’y(t)cos(cD j t)dt (10) 


Each sine wave must have an integral number of cycles over the 
entire run length (no partial waves), and no sine wave can have a 
frequency that is an integral multiple of another frequency. The 
relative amplitudes are selected so that the resulting input 
represents an input disturbance which occurs naturally in the task. 
Another benefit of using fourier coefficients is the Fast Fourier 
Transform (FFT). The FFT takes advantage of the periodic properties 
of sinusoids to reduce the computation time dramatically, but 
requires the number of data points to be an integer power of two. 

The least square error (LSE) method is the simplest approach 
mathematically, but computationally requires a large amount of 
storage space due to the matrix manipulations involved. In order to 



perform a least squares identification on a single-loop pilot-vehicle 
system an appropriate model must be chosen. Model selection varies 
from the simple crossover model, 



Ke'^ 

S 


(11) 


to higher order models like the precision model for single loop 
systems, (see ref. 2 pg 29) 



Once a model has been chosen, the input and model output error are 
minimized using a least squares technique. After the error has been 
minimized the coefficients are obtained, and the closeness of fit is 
determined. 

2.2 Least squares and sum of sines - Hess/Mnich. 

The Hess/Mnich research consisted of the identification of 
pilot dynamics from inflight tracking data using least squares and 
Fourier transform analysis. The NASA Ames Dryden Flight Research 
Facility with NASA Langley provided data from two flight tests for 
evaluation of pilot characteristcs. The task used for generating the 
data was an F-14 aircraft pursuing a T-38 target aircraft in both 
level flight and in a "3-G H wind up turn at a mach number of 0.55 at 
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an altitude of 10,000 feet with a separation distance of 800 feet. 
The F-14 pilot was using a gunsight reticle on a head-up display. 

The task of the F-14 pilot was to keep the reticle centered on the T- 
38 aircraft throughout the run. In addition to normal disturbances, 
the reticle in the F-14 was driven using a sum of sines as input so 
that the FFT results can be compared with the least squares results, 
(see ref. 6) 

The least squares technique used by Hess/Mnich for analysis 
is implemented in a software package called Nonintrusive 
Parameter Identification Program (NIPIP) (ref. 7). This program 
uses a general model with undetermined coefficients and determines 
the coefficients by comparing the data to the output of the model 
using a multiple linear regression technique (running least squares 
estimation). This program is capable of identifying linear and 
nonlinear relations between input and output as long as the 
relationships are linear with respect to the unknown coefficients. 
NIPIP uses a time frame length, the period over which the 
identification is to be performed, that can be specified as any part 
of the time history. NIPIP also has the option of sliding the frame 
along through the time history, removing old data as new data is 
entered, which yields a moving average through the time history. 

The sliding time window was not required for the Hess/Mnich 
analysis. 

The mathematical basis for the NIPIP program is a running 
least squares estimation technique. The coefficients of a prescribed 
difference equation approximating the relation ship between the 
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input and output of a linear system are estimated. For example, 
consider the following difference equation 

Y k = a 1 Y k-1 +b 1 x k + v k O 3 ) 

where Y k , X k , and V k are the system output, input, and modeling error at the k th 

sampling instant, respectively. Now considering a set of N measurements of the 
variables Y k and X k , one can write 

Y = HC' + V (14) 

where 



Now C' is found by minimizing the sum of the squares of V , where, 


V = Y - H C* 


(16) 
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and the sum of the squares is 

J = (Y - HC’) T (Y - HC’) (17) 

Minimizing a scalar J with respect to a vector C' requires 


aj 

9C'- 


0 


and 



(18) 


(19) 


applying Eqs. 18 and 19 to Eq. 17 yields 

H t HC' = H t Y (20) 

Solving for C' yields 

C’ = (H t H)- 1 H t Y (21) 

It can be shown that 
N 

hTh - X F k F k (22) 

k=1 


and 



10 


N 



(23) 


Then 


C’ = 



(24) 


where N is the number of data points, (see ref. 7) 

In addition to using NIPIP and the FFT to analize flight data 
Hess/Mnich also used a model to generate simulated data to test the 
two methods. Since the NIPIP program cannot identify time delays 
exactly, they chose their model with a second order denominator and 
a time delay that represents simplified human neuromuscular 
dynamics. The exact form of the mathetical model chosen was, 


2e* ,15s 

p= s 2 

<To> +04s+1 


(25) 


However, by changing the order of the model .time delays with 
integer multiples of the sampling period can be assumed and 
identified using a least squares technique. Then the quality of fit 
can be compared. This is the procedure followed in the Hess/Reedy 
research implemented in the Reedy subroutines. 
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Hess/Mnich produced good results with the simulated data. 
However due to a problem with the sum of sines input to the head up 
display the flight test results were not as well behaved as had been 
expected. Through averaging they were able to save the data and the 
results indicated that the crossover model fit the data in the area of 
crossover, (see ref. 6 fig. 9-12) 

2.3 Desciption of the Reedy identification - similarity 
to NIPIP 

The Reedy program consists of a least squares method applied 
to the data, either simulated or measured, where the transfer 
function is determined and converted from the z-domain into the w'- 
domain via the bilinear transform; 

1 +(T/2)w* 

2 " 1-(T/2)w* < 26 ) 

where T is the sampling rate. Since the coefficients of the discrete 
transfer function are nearly impossible to interpret in the discrete 
time domain, the Bode plots are used to convert to the frequency 
domain where analysis can be readily accomplished. The Bode plot is 
made using the w'-plane transfer function as an approximation of 
the S-plane. This approximation is valid when , (see ref. 8 pg 196) 

2 “s 

CO«t = 

1 n 


(28) 
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The previously described process is performed by two 
computer programs implemented as macros on CTRL-C, a computer- 
aided control system design package (ref. 9). The first macro is an 
identification of the coefficients of a difference equation 
representing the pilot model. The coefficients are determined by 
using a CTRL-C least squares method similar to the one described in 
part B of the background section. The identification macros are 
created by selecting a model from table 1; for example entry 3, 


where V is the output. X is the input, and b, and a, are the 
coefficients to be identified. Next the difference equation is found; 

Y(z) (1 -a, z' 1 ) - X(z) (b,z ’) (29 ) 


or, in the discrete time domain 


Y k-a,Y k ., -b,X k ., 


(30) 


then 


Y k " a i Y k-i + b i x k-i (31) 

The set of N measurements yield eqn (14), with the desired 

coefficients obtained from eqn (24). After the parameters have been 
determined the model is simulated to obtain the model output YV 

K 
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which is then compared with the original output Y k to obtain the 
quality of fit; 


r2 1 I>k - Y'k > 2 

IV 


(32) 


an R 2 value of unity indicates an exact fit. This procedure is 
Similar to the least squares portion of the NIPIP program without 
the sliding window. Although NIPIP would have performed the task, 
the least squares routines were written so that the student 
investigator would have a better understanding of the identification 
process and not simply be executing a "canned" program. 

The second macro has several characteristics that must be 
taken into account when operated independently of the first macro. 
First, the numerator and denominater must be defined as Num and 
Den respectively before running the transformation program. Num 
and Den must be the same size, they must be row vectors, and then 
coefficients must be in descending powers of 2. If they are not the 
same order or contain zero coefficients, zeros must be added 
appropriately inorder to achieve this constraint. The transformation 
macro operation is very simple when used in conjuncton with one of 
the compatible least squares identification macros. The only 
information necessary is the sampling rate. The program will pause 
to display the W -plane transfer function, and after depressing the 
return key the magnitude and phase Bode plots will be constructed. 
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Both macros are completely self contained programs and can be ran 
independently when necessary. Copies of these programs are in 
appendix III. 

2.4 Example cases; Noise effects, Input Bandwidth 
effects, Delay Estimation techniques. 

The example cases are simulations performed on the Advanced 
Continuous Simulation Language (ACSL) (ref. 10) to generate data for 
the identification process. Many examples were run to build an 
understanding of the identification process and gain experience with 
the procedures. The ACSL simulation programs were written by 
Ronald A. Hess, Professor of Mechanical Engineering University 
California at Davis. Several examples of the simulation programs 
are listed in appendix III. 

A total of six examples were run each varying in complexity. 
Five test cases were run using the same system transfer function, 
and one higher order multi-loop example was used. This transfer 
function 


1 0Oe'^ 8 
s[s 2 +4s+100] 


(33) 


was chosen because it represents a second order system with an 

integration and a time delay and is typical of pilot-vehicle dynamics 
(Y p V c ) in single-loop tasks. The simulations were run with a time 

step of 0.05 seconds for 102.4 seconds yielding 2048 data points 
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(2 11 ), thereby meeting the "power of 2" requirement for FFT analysis. 
The FFT was used as a comparison to the least squares in some of 
the examples. The first five test cases run were, an open loop 
system with To = 0-0 seconds (figure 1), a closed single loop system 
with ib-O.O seconds (figure 2), a closed single loop system with 
lo^O.O seconds and injected noise (figure 3), a closed single loop 
system with To=0.3 seconds and injected noise, and a closed single 
loop system with ib*0.0 seconds, injected noise and a bias error 
(figure 4). In each of these test cases the C(s)/E(s) transfer 
function was identified. The last example was a multi-loop hovering 
helicopter with a realistic pilot model with noise and a bias error 
(figure 5). The Q/Xe transfer function was identified in this case. 
The Q/Xe transfer function is a "composite", and it is similar to a 
transfer function to be measured in the CH-47B simulation to be 
discussed later. From these six examples many comparisons can be 
made. The effects of noise, bias, delay, cutoff frequency, and open 
loop versus closed loop dynamics on the quality of the 
identifications will be discussed presently. 

An injected noise signal with a root-mean-square (RMS) value 
of 0.1 times the input RMS does not effect the identification to any 
appreciable extent. This is indicated in figure 6 and table 2. 

Increasing the RMS value to 1.0 times the input RMS, however, 
severely compromises the least- squares identification. 

In table 2 and in subsequent tables listing the model used for 
the identification of the transfer function, the model column 
contains a code. The first two digits represent the entry position in 



table 1 ,the second two digits represent the time delay in integer 
multiples of the sampling rate, and the last digit indicates if the 
routine contains a bias identification. For example, Z8D0B 
represents entry 8 with zero delay and a bias identification. 

The effect of bias appears in figure 7 and table 3. Although 
not as dramatic as the noise effects, the bias error reduces the 
ability to identify the system correctly. In table 3 the last three 
rows correspond to runs in which a unity constant bias term was 
added to the simulation (fig. 4), and when an identification routine 
with a bias identification is used a definite improvement results. 

The effect of delay is apparent in figure 8 and tables 4 and 5. 
As expected, increasing the delay decreases the correlation 
coefficient. It is important that the delay be correctly identified 
for best results. Table 5 shows that when the delay in the 
identification matches the delay in the system the correlation 
coefficient is maximized. 

The effect of input cutoff frequency is shown in figure 9, 10 
and 1 1 and in tables 2, 3, 4, and 6. Cutoff frequency appears to have 
no effect on the open loop and closed loop with zero noise, zero bias 
and zero delay examples. However, in ail other test cases the 
increase in cutoff frequency reduced the correlation coefficient. 

This effect is most noticeable in the delay, noise, zero bias 
example, and can be seen in the last three rows of table 4. The sum 
of sines cutoff frequencies and their corresponding magnitudes are 
listed in table 7. 
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Identification of the open-loop transfer function in closed- 
loop versus open-loop operation was also undertaken and the results 
appear in figure 12 and 13 and in table 6. Two different inputs were 
utilized; the sum of sines and a unit step. The sum of sines, open- 
loop identification varies slightly at the high and low ends of the 
frequency range, but the deviation is insignificant. Identification of 
the open-loop transfer function in closed-loop versus open-loop 
operation with zero noise using both FFT and LSE techniques was 
undertaken. The medium cutoff sum of sines input was used. 

The results appear in figures 14 and 15. 

The last example was the multi-loop hovering helicopter with 
a realistic pilot model. Comparisons were made between the FFT, 
the least squares without bias, and the least squares with bias. The 
results show that the least squares with bias corresponds very well 
with the FFT, and yields a high correlation coefficient. See figures 
16, 17 and table 8. The value of bias calculated by the identification 
subroutines cannot be compared with the bias in the single loop 
examples, due to the location in the model were the signals are 
measured. 

3.0 A multi-axis manned simulation task 

Data obtained from the CH-47B variable-stability helicopter at 
the NASA Ames Research Center was analyzed using the LSE and the 
FFT identification procedures outlined previously. The transfer 
functions identified are presented in bode form to facilitate 
comparison. The experimental tracking task was performed while 
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the helicopter was in the attitude command/attitude hold dynamic 
mode (ref. 11). During the precision tracking task the pilot 
attempts to maintains a hovering position above a pad symbol, while 
the pad symbol is driven by a forcing function. The forcing function 
is a random appearing sum of sinusoids (tables 9 and 10). The 
sampling interval for the simulation was 0.05 seconds. A run length 
of 102.4 seconds yielding 2048 (2 11 )data points was used in order 
to meet the "power of 2" requirement of the FFT analysis technique. 

The subjects used in this experiment were an engineer and a 
test pilot. Each of the two subjects performed five runs on three 
different configurations of a tracking task. In the first 
configuration the pilot controled only the longitudinal motion of the 
helicopter with the aid of the complete display (figure 18 ). In the 
second configuration the subject controlled both the lateral and 
longitudinal motion of the helicopter with the aid of the complete 
display. In the final configuration the subject controlled both the 
lateral and longitudinal motion of the helicopter, but without the 
acceleration symbol on the display. 

Three transfer functions were analyzed using the different 
configurations; these were: X/Xpad, Xball/Xpad, and Xball/Xball error 
(see figure 19). The X/Xpad transfer function is defined between the 
longitudinal position of the helicopter and the longitudinal position 
of the hover pad, measured in feet from a fixed point on the earth. 

The Xball/Xpad "composite" transfer function is defined between the 
position of the acceleration symbol relative to the center of the 
display and the position of the hover pad measured from a fixed 



point on the earth, both measured in display units. The Xball/Xball 
error transfer function is defined between Xball and the 
longitudinal error between the displayed hover pad and the 
acceleration ball, again, both in screen units. The X/Xpad transfer 
function was analyzed using all three configurations, while the 
composite and the Xball/Xball error transfer functions were 
analyzed using only the first two configurations. 

For the X/Xpad transfer function identifications a good 
correlation exists between the LSE and the FFT when the 
acceleration symbol is present (see figures 20 - 23 ), but when the 
acceleration symbol is removed the comparison becomes poor (see 
figures 24 and 25 ). The correlation coefficients are lower for LSE 
identifications without the acceleration symbol (see tables 11-16 
), and a lightly-damped mode appears in the bode plots. The mode, 
evident in both the FFT and the LSE plots, indicates a decrease in 
closed-loop system stability. Another comparison demonstrating 
the utility of the acceleration symbol appears in figures 26 and 27 
These two figures are plots of the actual output X, the simulated 
output X’ for subject 1 and the command signal Xpad. The simulated 
output comes from the LSE program just before calculating the 
correlation coefficient. After the parameters or coefficients of the 
transfer function are identified, the program simulates the 
identified transfer function and compares the simulated output to 
the actual output. It is this simulated data that is plotted with the 
actual data versus time in figures 26 and 27. The correlation 
coefficients for figure 26 and 27 are 0.9883 and 0.7997 respectively 
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(see tables 13 and 15 ). From this analysis it is obvious that the use 
of the acceleration symbol greatly increases the pilots ability to 
maintain a hovering position over a moving object. The poor 
performance evident in figure 26 even with the acceleration ball is a 
result of the challenging nature of the sum of sines input. 

Figures 28-31 show the LSE and FFT measurements for the 
Xball/Xball error transfer function, for each subject and 
configuration. Tables 21-24 show the pertinent parameters for 
these measurements. As the tables indicate, only the FFT 
measurements were reliable for this transfer function. The LSE 
technique yielded either unstable transfer functions (unbounded 
values) or very low values. This poor identification performance 
with the LSE technique may be do to the effects of noise injection by 
the subjects (remnant). Note that large noise injection did adversly 
effect LSE identification performance in the example case in Section 
2.4. 

The FFT results of figures 28-31 and tables 21-24 were quite 
acceptable. The data can be interpreted in terms of the crossover 
model of Eq. 11, with crossover frequencies on the order of 2.0 
rad/sec and time delays of approximately 0.2 secs. 

The composite transfer function Xball/Xpad yields important 
information about the assumed pilot control structure shown in 
figure 19. This loop structure assumes single-loop compensatory 
behavior on the part of the pilot, i.e. that Xe, itself, is not used by 
the pilot, only what has been called Xball error. Thus, although 
figure 19 shows two loops being closed, only the inner loop is 
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assumed to be closed by the pilot. Now, in order to improve tracking 
performance , the pilot might adopt what would be interpreted as 
pursuit behavior in the single loop manual control structure of 
figure 19 (ref. 2), or multi-loop behavior in terms of the multi-loop 
manual control structure of figure 32. In other words Xe might be 
used by the pilot and be subject to compensation. The resulting 
compensated signal would then be compared with Xball and the 
difference be used to close the inner control loop. This is the multi- 
loop manual control structure shown in figure 32. Note that the 
inner loop error signal is now not Xball error, but some internally 
generated error based upon the difference between the compensated 
Xe (called Xec in figure 32 ) and Xball. 

Figures 33-36 show measured composite transfer functions 
Xball/Xpad for the two subjects for longitudinal tracking alone and 
simultaneous longitudinal and lateral tracking. Now the transfer 
function Xball/Xpad can be writen as a product or composite of two 
other transfer functions as, 

x ball x ball x e 

x pad x e x pad 

Now, the second of these, Xe/Xpad, is the error-to-input transfer function for the 
outer loop. Figures 20-23, show the closed loop X/Xpad transfer functions 
obtained in this study. Defining bandwidth as that frequency where the phase 
goes through -90 degrees these figures indicate a closed loop bandwidth of 
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around 0.30 rad/sec. This means that the Xe/Xpad transfer function will be very 
close to unity for all frequencies much beyond 0.30 rad/sec. Thus 


x ball ^ x ball 
x pad x e 


(35) 


for to > 0.30 rad/sec. But the transfer function can be obtained as 


x ball x bal/ x ball error 

x e 1 + X| 3a ||/X^ a || error 


(36) 


Recall that the Xball/Xball error transfer function has already been obtained, at 
least in terms of FFT measurements (see figures 28-31). Taking the FFT 
measurements of figure 28 as a representative sample, an acceptable fit to the 
data was obtained in the form of a rational transfer function. Now, forming 
Xball/Xe using this fit, and approximating Xe/Xpad as 

X e s 
x pad s + 0-3 

one obtains figure 37 from the product on the right hand side of Eq. 34. If no 
compensation of Xe is occuring, then figure 34 should resemble figures 33-36 
in the frequency range to > 0.03 rad/sec. Looking at the amplitude ratio, 
this is not the case for frequencies above around 2.0 rad/sec. The 
measurement of figures 33-36 all indicate that the amplitude ratios 
are relatively flat and greater than unity in value, whereas the 
amplitude ratio of figure 37 has begun to fall off at around 2 



rad/sec. The LSE result with the highest R 2 value (run 4, table 20) is 
shown for comparison. Thus, some form of pilot compensation as 
suggested in figure 32, is probably occuring in the outer loop to 
cause this discrepancy. This does not imply that the FFT 
measurements assuming the loop closure structure of figure 19 are 
incorrect, rather they reflect the effective compensatory behavior 
of the pilot. 

Thus the FFT and least-squares measurements of the 
composite transfer function have led to the discovery of pursuit 
behavior in terms of a single-loop manual control structure as 
shown in figure 19, or, equivalently, multi-loop behavior, in terms 
of the control structure of figure 32. While the evidence for this 
behavior has been obtained indirectly, the data supporting it has 
been quite consistent. As the pertinent figures and tables for the 
Xball/Xpad transfer functions indicate, twice the standard deviation 
of the FFT data is typically less than a symbol width in magnitude, 
and the R 2 values for the least squares data are typically greater than 0.95. 
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4.0 Conclusions 

1) After preliminary investigation with synthetic data Fast 
Fourier Transform (FFT) and Least Square Error (LSE) estimation 
techniques were applied to the identification of pilot-vehicle 
dynamics in a realistic flight simulation task. 

2.) With the exception of the identification of the Xball to 
Xball error transfer function, comparisons between the FFT and LSE 
techniques were, in general, good. No acceptable LSE identification 
of the aforementioned transfer function was found. It was thought 

this poor LSE performance might be attributed to human noise 
injection in the inner control loop. 

3. ) The FFT identification of the Xball to Xball error transfer 

function could be described in terms of the well-known crossover 
model of the human pilot. 

4. ) The Identification of a composite transfer function 

yielded some indirect evidence of pursuit tracking behavior on the 
part of the test subjects. 
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Tables 
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TABLE 1 


z-T ransform 


b,*’ 1 


b,Tz-' 

1 - z' 1 

b i z '’ 

b^'UbgZ* 2 

1 -a 1 z* 1 -a 2 z* 2 


No. of 
Unknowns 


Physical System 

K (pure gain) 

K 

7* (integrator) 


Ka 

s+a (ist-order lag) 


1 tu 2 4 

l-a^^-agZ' 2 ^ 

4 

K(s+a) 

roj 

b^+b-z* 2 
l-a.z' 1 (Z ' P) 

3 

K(s+a) 

(s+b) 

b 1 z‘ 1 +b 2 z* 2 
1 -a 1 z' 1 -a 2 z‘ 2 

4 

K(s+a) 

s(s+b) 

b 1 2 ' 1 +b 2 z' 2 +b a z' 3 

1 - a <*-1 - --2 _ _.3 

6 

K(s+a 


b 1 z' 1 +b 2 z- 2 +b 3 z- 3 +b 4 z- 4 

~ l-a^-’-agZ^-agZ ' 3 


K Kn*> n ] 



TABLE 2 - NOISE AND CUTOFF FREQUENCY COMPARISON 


ACTUAL 

INPUT R 2 NOISE DELAY 


ZSINES 




MED CUT 

0.9997 

0.0 

0.0 

XSINES 

LOW CUT 

0.9940 

0.1 

0.0 

ZSINES 

MED CUT 

0.9809 

0.1 

0.0 

ZSINES 

HIGH CUT 0.9552 

0.1 

0.0 

ZSINES 
MED CUT 

0.1114 

1.0 

0.0 


IDENTIFIED LOOP 
BIAS STRUCTURE 

MODEL 

N.l. 

CLOSED 

28D0 

N.l. 

CLOSED 

28D2 

N.l. 

CLOSED 

28D2 

N.l. 

CLOSED 

28D2 

N.l. 

CLOSED 

28D4 


* N.l. - NOT IDENTIFIED 



TABLE 3 - BIAS COMPARISON 


INPUT* R 2 


ACTUAL 

IDENTIFIED LOOP 

NOISE 

DELAY 

BIAS 

STRUCTURE 

XSINES 

LOW CUT 0.9802 

XSINES 

0.1 

0.3 

N.l. 

CLOSED 

MED CUT 0.8916 

XSINES 

0.1 

0.3 

N.l. 

CLOSED 

HIGH CUT 0.7276 
XSINES 

0.1 

0.3 

N.l. 

CLOSED 

LOW CUT 0.9969 
XSINES 

0.1 

0.3 

0.9961 

CLOSED 

MED CUT 0.9904 
XSINES 

0.1 

0.3 

0.9911 

CLOSED 

HIGH CUT 0.9839 

0.1 

0.3 

0.9954 

CLOSED 


MODEL 

Z8D4 

Z8D4 

Z8D5 

Z8D3B 

Z8D3B 

Z8D4B 


* N.l. - NOT IDENTIFIED 
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TABLE 4 - DELAY COMPARISON 


ACTUAL IDENTIFIED LOOP 

INPUT R 2 NOISE DELAY BIAS STRUCTURE 


ZSINES 
LOW CUT 

0.9940 

0.1 

0.0 

N.L 

CLOSED 

ZSINES 
MED CUT 

0.9809 

0.1 

0.0 

N.L 

CLOSED 

XSINES 

HIGH CUT 

0.9552 

0.1 

0.0 

N.L 

CLOSED 

XSINES 

LOW CUT 

0.9802 

0.1 

0.3 

N.L 

CLOSED 

ZSINES 
MED CUT 

0.8916 

0.1 

0.3 

N.L 

CLOSED 

ZSINES 
HIGH CUT 

0.7276 

0.1 

0.3 

N.L 

CLOSED 


MODEL 

28D2 

28D2 

28D2 

28D4 

28D4 


* N.l. - NOT IDENTIFIED 


28D5 



TABLE 5 - IDENTIFIED DELAY COMPARISON 
ACTUAL IDENTIFIED LOOP 

INPUT R 2 NOISE DELAY BIAS STRUCTURE MODEL 


XSINES 

MED CUT 

0.9338 

0.1 

0.0 

XSINES 

MED CUT 

0.9540 

0.1 

0.0 

SSINES 
MED CUT 

0.9809 

0.1 

0.0 

XSINES 

MED CUT 

0.9555 

0.1 

0.0 

ZSINES 
MED CUT 

UNSTABLE 

0.1 

0.0 

XSINES 

MED CUT 

UNSTABLE 

0.1 

0.0 

XSINES 

MED CUT 

UNSTABLE 

0.1 

0.0 


N.l. 

CLOSED 

28D0 

N.l. 

CLOSED 

Z8D1 

N.l. 

CLOSED 

Z8D2 

N.l. 

CLOSED 

Z8D3 

N.l. 

CLOSED 

Z8D4 

N.l. 

CLOSED 

Z8D5 

N.l. 

CLOSED 

Z8D6 


* N.l. - NOT IDENTIFIED 
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TABLE 6 - OPEN VS. CLOSED LOOP 


ACTUAL IDENTIFIED LOOP 


INPUT R 2 

NOISE 

DELAY 

BIAS 

STRUCTURE 

MODEL 

XSINES 

LOW CUT 1 .0 

0.0 

0.0 

N.l. 

OPEN 

Z8D0 

Esines 

MED CUT 1.0 

0.0 

0.0 

N.l. 

OPEN 

Z8D0 

Esines 

HIGH CUT 1.0 

0.0 

0.0 

N.l. 

OPEN 

Z8D0 

XSINES 

LOW CUT 0.9998 

0.0 

0.0 

N.l. 

CLOSED 

Z8D0 

XSINES 

MED CUT 0.9997 

0.0 

0.0 

N.l. 

CLOSED 

Z8D0 

XSINES 

HIGH CUT 0.9998 

0.0 

0.0 

N.l. 

CLOSED 

Z8D0 

STEP 1 .0 

0.0 

0.0 

N.l. 

OPEN 

Z8D0 

STEP 0.9999 

0.0 

0.0 

N.l. 

CLOSED 

Z8D0 

* N.l. - NOT IDENTIFIED 







TABLE 7 - XSINES 


12 

* A oX A i sin ^ a) i t ) 

i=i 



LOW CUTOFF 
A q =1/1 .546 

MED CUTOFF 
A 0 «1/3.03 

HIGH CUTOFF 
A 0 = 1/4.02 

HELICOPTER 
A n = 1 

i 

COj (rad/sec) Aj 

A i 


A , 

1 

0.1841 

1 

1 

1 

17.6 

2 

0.3068 

1 

1 

1 

17.6 

3 

0.4909 

1 

1 

1 

17.6 

4 

0.7977 

0.1 

1 

1 

17.6 

5 

1.1660 

0.1 

1 

1 

1.76 

6 

1 .7790 

0.1 

1 

1 

1.76 

7 

2.8230 

0.1 

0.1 

1 

1.76 

8 

4.6630 

0.1 

0.1 

1 

0.88 

9 

6.9330 

0.1 

0.1 

0.1 

0.88 

10 

8.9580 

0.1 

0.1 

0.1 

0 

11 

12.0880 

0.1 

0.1 

0.1 

0 

12 

17.9780 

0.1 

0.1 

0.1 

0 



TABLE 8 - MULTI-LOOP HELICOPTER COPARISON 


INPUT 

R 2 

NOISE 

XSINES 

0.9644 

0.0 

Zsines 

0.9664 

0.0 

ZSINES 

0.9573 

0.1 

XSINES 

0.9578 

0.1 


ACTUAL IDENTIFIED LOOP 

DELAY BIAS STRUCTURE MODEL 


0.0 

0.0 

0.0 

0.0 


N.l. 

-14.01 

N.l. 

-13.95 


MULTI 

MULTI 

MULTI 

MULTI 


Z9D0 

Z9D0B 

Z9D0 

Z9D0B 


* N.l. - NOT IDENTIFIED 
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47B HELICOPTER S F ° R ™ E LONGITUDINA1 - DIRECTION OF THE 

12 

r ( t ) ■ A oX A i sin ^ OJ i t ) 

i-1 

ct)j (rad/sec) AMPLITUDE NO. OF CYCLES 


1 

0.1841 

17.6 

3 

2 

0.3068 

17.6 

5 

3 

0.4909 

17.6 

8 

4 

0.7977 

17.6 

13 

5 

1.1660 

1.76 

19 

6 

1.7790 

1.76 

29 

7 

2.8230 

1.76 

46 

8 

4.6630 

0.88 

76 


9 6.9330 


0.88 


113 
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TABLE 10 - XSINES F0R THE 
HELICOPTER. 


LATERAL DIRECTION OF THE CH-47B 


12 

r(t > " A oX A i S ' n(a) l t ) 

i-1 


1 

toj (rad/sec) 

AMPLITUDE 

NO. OF CYCLES 

1 

0.2454 

17.6 

4 

2 

0.4295 

17.6 

7 

3 

0.6750 

17.6 

11 

4 

0.9204 

17.6 

15 

5 

1.4110 

1.76 

23 

6 

2.2700 

1.76 

37 

7 

3.7430 

1.76 

61 

8 

5.7060 

0.88 

93 

9 

7.7930 

0.88 

127 



38 




TABLE 1 1 
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TABLE 12 
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TABLE 14 





in 

LU 
I 

CD 

< 

h- 













































































co 


ca 

-Q 


43 


c 



1 

£00 

100 

-522 

O 

r- 

g 

co 

L_ 

d> 

CD 

O 

O 

si 

910 

900 

0.02 

-461 

r- 

Cu 

O 

c 

CD 

c 


0.12 

oo 

O 

d 

O 

o 

O 

1 

o 

o 

CO 

w. 

4-^ 

75 

c 

in 


990 

0.41 

930 

-175 

oo 

ON 

FFT results for the X/Xpad transfer function, subject 2, lateral and longitu< 

to VO 

3 Zi 

1.07 

190 

0.35 

ON 

1 

r^ 

B 

K 

0.34 

0.26 

O 

cs 

o 

o 

1 


l 

E 

K 

0.62 

0.41 

0.27 1 

vO 

7 

oo 

oo 

q O 
^ c<i 

o 

1.42 

V-l 

OO 

d 

0.52 

r- 

1 

r- 

04 

©i 

0.1841 

1.50 

0.84 

0.47 

-68 

CO 

I 

0 

+ 

p 

P 

0 

1 

p 



1 

MAGNITUDE 

PHASE 



COEFFICIENTS OF THE IDENTIFIED TRANSFER FUNCTION IN 
THE Z-DOMAIN IN DESCENDING ORDER NUMERATOR FIRST 

2 . 9957 , - 2 . 9936 , 0 . 9979 ; 0 . 0012 , - 0 . 0027 , 0.0015 


2 . 9888 , - 2 . 9815 , 0 . 9927 ; 0 . 0003 , - 0 . 0007 . 0.0004 

2 . 9912 , - 2 . 9857 , 0 . 9945 ; - 0 . 0001 , 0 . 0003 , - 0.0001 

2 . 9889 , - 2 . 9826 , 0 . 9937 ; 0 . 0006 , - 0 . 0012 , 0.0006 

ENTIFIED 

DELAY 

CM 

1 

N 


‘k 

Ssi 


o 


■ 




o 


■ 




LU 

CM 


o 

ca 


LL CO 

Tj " 


CO 


T— 

F < 

in 

■ 

05 

in 

CJ 

Z CD 

C ". 


y— 

CO 


UJ 

Q 

d 


d 

o 

LO 


H 


l 

1 

i 


00 

LU 

_J 

CO 

O) 

CD 

CO 

CNJ 

05 

< 

05 

CO 

o 

cc 

00 

h - 

T— 

y— 

o 

O 

CO 


* t 

00 


• 

z 

• 

• 



o 


O 

O 

o 

Nna 


CNJ 



in 


c 

o 

75 

0) 

<D 

O 

o 

co 

o 

c 


O) 

I* 

o 

CO 


CO 

E 

T3 

3 

■*-* 

O) 

c 

o 


"O 

c 

(0 

75 

L. 

3 

co 


CNJ 

O 

<D 

lo 

3 

</) 


c 

o 

4^ 

o 

c 

3 


3 

7/5 

c 


CO 


“D 

co 

CL 

X 

*: 


0 


</> 

<4—* 

3 

(/) 

CD 


CD 

HI 

_l 

CD 

< 

I- 


LU 

C/5 




















































44 



ransfer function, subject 1, longitudinal tracking only. 
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ransfer function, subjoct 2, latoral and longitudinal tracking. 
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FIGURE 2 - CLOSED LOOP 



FIGURE 3 - CLOSED LOOP WITH NOISE N(s) 



FIGURE 4 - CLOSED LOOP WITH NOISE N(s) AND BIAS B(s) 



FIGURE 5 - MULTI-LOOP HOVERING HELICOPTER MODEL 
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FIGURE 8 - DELAY COMPARISON, MEDIUM CUTOFF FREQUENCY. 
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FIGURE 16 - MULTI-LOOP FFT VS. LEAST SQUARES WITH BIAS, AND ZERO NOISE 
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FIGURE 17 - MULTI-LOOP FFT VS. LEAST SQUARES WITH NOISE AND BIAS 
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WARMING 


FIGURE - 18 Hover Display Symbology. 
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FIGURE - 19 Hypothesized pilot loop closures with acceleration ball display 
symbology. 
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FIGURE 30 - Xball/Xball error TF ID for Subject 1, Lateral and longitudinal trackina 
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31 - Xball/Xball error TF ID for Subject 2, Lateral and longitudinal tracking. 
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FIGURE 35 - Xball/Xpad TF ID for Subject 1, Lateral and longitudinal tracking. 
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8.0 Appendix III 


Computer Programs 



** ACSL multi-loop helicopter simulation 

PROGRAM HEU1 

CONSTANT TFIN-1 10. 

CONSTANT TDX-0.15 
CONSTANT KEX1-1. 

CONSTANT QX2-1..14.14, 100. 
CONSTANT K-98.03 
CONSTANT K1-1. 

CONSTANT K2-10. 

CONSTANT T1-2.5 
CONSTANT T2-.369 
CONSTANT XU-0.05 
CONSTANT G-32.2 
CONSTANT MU-.0311 
CONSTANT MTH-8. 

CONSTANT MDE-.585 
CONSTANT MQ-5. 

CONSTANT TL-3.0 
CONSTANT KY-.0071 
CONSTANT ta-0.2 
CONSTANT kn-1.0 
CONSTANT m-0.0 
CONSTANT s-0.02710 
CONSTANT mn-0.00121285 

ARRAY QX2(3) 

CINTERVAL CINT-0.05 

INITIAL 

UNIFI(1717) 

END $ "OF INITIAL" 

DYNAMIC 

DERIVATIVE 



DX=TRAN(0,2,100.,QX2,U1) 

UM 1 =(.01 *T 1 )*LEDLAG(1 00., T1 ,DX,0.) 

"UMX1=K2*UM1" 

"UMX1 =(K2/T2)*LEDLAG(T2,.01,UM1,0.)" 

UMX1 =(K2*T2)*LEDLAG(0.,T2,UM1 .0.) 

U1 =UCX1 -UM1 -UMX1 

UCX1 =KEX1 *DELAV(X1 EN,0.,TDX,300) 

xle =thc-th 

XI EN»THC-TH+en+0.1 
en=kn*(ou(ta,m,s)-mn) 

XC=1 7.6*(SIN(.1 841 *T)+SIN(.3068*T)+SIN(.4909*T)+SIN( 7977‘T) 
0*1 *(SIN(1 .1 66*T)+SIN(1 .779*T)+SIN(2.823*T))+ 
0.05*(SIN(4.663*T)+SIN(6.934*T))) 

XE=XC-X 

THC— KY*LEDLAG(TL,.01 ,XE,0.) 

UD=XU*U-G*TH 

THD=Q 

QD=MU*U+MTH*TH+MQ*Q+MDEVDX*K) 

XD=U 

TH=INTEG(THD,0.) 

Q=INTEG(QD,0.) 

U=INTEG(UD,0.) 

X=INTEG(XD,0.) 

TP-1./TFIN 

VX1 E»TP*INTEG(X1 E**2.,0.) 

VTHC=TP*INTEG(THC**2.,0.) 

VXC=TP*INTEG(XC*"2.,0.) 

VX=TP*INTEG(X**2.,0.) 

VDX=TP*INTEG(DX**2.,0.) 

nbar»tp*integ(en, 0 .) 

END $ M OF DERIVATIVE" 

TERMT(T.GE. TFIN) 

END $ "OF DYNAMIC" 

END $ "OF PROGRAM" 



** ACSL closed-loop simulation, sum of sines input with noise, ** 

** delay , and bias. 

program ndbsin 
array d(4) 

constant d-1.,4.,100.,0., ta-0.2, kn-1.0, mn-0.00236539, m»0., s-0.03215 

constant wl -0.1 841, w2«0.3068, W3-0.4909, W4-0.7977, w5-1.166 

constant W6-1.779, W7-2.823, W8-4.663, W9-6.933, wlO-8.958 

constant wl 1 -1 2.088, wl 2-1 7.978, td-0.3 

tfin-102.4 

cinterval cint-0.05 

INITIAL 

UNIFI(3333) 

END $"OF INITIAL" 

dynamic 

derivative 

ed=delay(e1 +n,0.,td,300) 

c=tran(0,3,100.,d,ed) 

el-r-c+1. 

e-el-1. 

n=kn*(ou(ta,m,s)-mn) 

r»(sin(wrt)+sin(w2*t)+sin(w3*t)+0.1*(sin(w4*t)+sin(w5*t)+... 
sin(w6*t)+sin(w7*t)+sin(w8*t)+sin(w9*t)+sin(w10*t)+... 
sin(w1 1 *t)+sin(w1 2*t)))/1 .546 
tf-1 ./tfin 

nbar-tf*integ(n,0.) 
end $ "of derivative" 
end $ "of dynamic" 
termt(t.ge.tfin) 
end $ "of program" 
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** ACSL closed-loop simulation, sum of sines input with noise, ** 
and delay . 


program ndsin 
array d(4) 

constant d«1.,4.,100.,0., ta-0.2, kn-1.0, mn-0.00303933, m«0.0, s-0.04131 

constant wl -0.1 841, W2-0.3068, W3-0.4909, W4-0.7977, W5-1.166 

constant W6-1.779, W7-2.823, W8-4.663, W9-6.933, wl 0-8.958 

constant wl 1 -1 2.088, wl 2-1 7.978, td-0.3 

tfin-102.4 

cinterval cint«0.05 

INITIAL 

UNIFI(3333) 

END $"OF INITIAL" 

dynamic 

derivative 

ed-delay(e+n,0.,td,300) 

c-tran(0,3,100.,d,ed) 

e=r-c 

n*kn*(ou(ta,m,s)-mn) 

r-(sin(w1*t)+sin(w2*t)+sin(w3*t)+sin(w4*t)+sin(w5*t)+... 
sin(w6*t)+0.1 *(sin(w7*t)+sin(w8*t)+sin(w9*t)+sin(w1 0 *t)+... 

sin(w1 1 *t)+sin(w1 2*t)))/3.03 
tf-1 ./tfin 

nbar-tf*integ(n,0.) 
end $ "of derivative" 
end $ "of dynamic" 
termt(t.ge.tfin) 
end $ "of program" 



** ACSL closed-loop simulation, sum of sines input with noise. ** 


program nsin 
array d(4) 

constant d-1 .,4., 1 00., 0., ta-0.1, kn-1.0, mn-0.05276090, m-0., s«1 0 

constant wl -0.1 841 , W2-0.3068, W3-0.4909, W4-0.7977, w5-1 1 66 

constant w6-1. 779, W7-2.823, W8-4.663, W9-6.933, w10-8 958 

constant wl 1 -1 2.088, wl 2-1 7.978 

tfin-102.4 

cinterval cint-0.05 

INITIAL 

UNIFI(3333) 

END $ H OF INITIAL" 

dynamic 

derivative 

c=tran(0,3,100.,d,e+n) 

e-r-c 

n-kn*(ou(ta,m,s)-mn) 

r=(sin(w1 *t)+sin(w2*t)+sin(w3*t)+sin(w4*t)+sin(w5*t)+... 
sin(w6*t)+0.1 *(sin(w7*t)+sin(w8*t)+sin(w9*t)+sin(w1 0*t)+. 
sin(w1 1‘t)+sin(w12*t)))/3.03 
tf=1 ./tfin 

nbar-tf*integ(n,0.) 
end $ "of derivative" 
end $ "of dynamic" 
termt(t.ge.tfin) 
end $ "of program" 



** ACSL closed-loop simulation with sum of sines input. 


** 


program sin 
array d(4) 

constant d«1 .,4., 100..0. 

CONSTANT w1-0.1841.w2-0.3068, w3- 0.4909. W4-0.7977. w5-1 166 

constant w6-1. 779, W7-2.823, W8-4.663, W9-6.933 

CONSTANT W1 0-8.958, W1 1 -1 2.088, W1 2-1 7 978 

tfin-102.4 

cinterval cint=0.05 

dynamic 

derivative 

c=tran(0,3,100.,d,e) 

e=r-c 

r»(sin(w1*t)+sin(w2*t)+sin(w3*t)+sin(w4*t)+sin(w5*t)+... 
sin(w6*t)+sin(w7*t)+sin(w8*t)+0.1*(sin(w9*t)+sin(w1 0 *t)+... 

sin(w1 1 *t)+sin(w1 2*t)))/4.02 
end $ "of derivative" 
end $ "of dynamic" 
termt(t.ge.tfin) 


end $ "of program" 



// *** 2 -domain to w*-domain transformation subroutine 
T-0.05; 

ZZ-SIZE(NU); 

N*2Z(1 ,2)-1 ; 

X -IT/21]; 

Y = [-T/2 1]; 

K-N+1; 

A=0*ONES(N,K); 

B=0*ONES(N,K); 

C-0;D-0;DD-0;CC-0;NUM-0;DEN-0; 

YY-1;M-1;XX-1; 

FOR 1-1 
FOR J-1 
L-J+K-l;... 

A(I,L)-XX(1,J);... 

B(I,L)-YY(1,J);... 

END;... 

XX-CONV(XX.X);... 

YY»CONV(YY,Y);... 

M-M+1 ;... 

END; 

Q-K; 

FOR P-1 :K;... 

C(P,:)-NU(1 1 P)*CONV(A(Q l :),B(P, 

D(P 1 :)-DE(1 1 P)*CONV(A(Q,:) l B(P i :));... 

Q-K-P;... 

END; 

CC-C';DD-D';S-2*N+1 ;Q2-K; 

FOR R-1 :K;... 

NUM(1 ,Q2)-SUM(CC(S, 

DEN(1,Q2)-SUM(DD(S, 

Q2-K-R;... 

S-S-1;... 

END; 

NUM.DEN 

n=NUM; 

q-DEN; 

PAGE 



[a,b,c,d]«tf2ss(n,q); 

v=logspace(-1,2); 

[mag,pha]«bode(a,b,c,d,1 ,v); 

W1ND0W('21 1') 

a»[0.1 ,0.01 ;1 00,1 0;33.3,3.33J; 

plot(a, 'scale') 

plot(v, mag, 'loglog', 'dotted') 

title(’magnitude',' ') 

WINDOW('212') 

aa-[0. 1 ,-300;1 00,-50:33.3,50]; 

plot(aa, 'scale') 

plot(v, pha.'logx', 'dotted') 

title('phase',' ') 
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LOAD FAST F0UR,ER TRANSFORM subroutine *** 

T=H(153:2200,1);’ 

E=H(1 53:2200,2); 

C=H(1 53:2200, 3); 

R=H(1 53:2200, 4); 

FFTCaFFT(C); 

FFTE=FFT(E); 

FFTR=FFT(R); 

W=[.1 841 ;.3068;.4909;.7977;1 .1 66-1 779- 
2.823 ;4.663;6.934J; 

FOR F*1 :9;.„ 

N*ROUND(W(F)*2048/(20*2*PI)+1)- 

MAG Em f AL<FFTC(N))) "^('MAGfFFTCW)) “2} 

PHAr/n^ 

FOr’e-1:9;,.. 

N-ROUND(W(Er2048/(20*2‘P|Ul V 
^IFREAL(FFTC(N)) < 0. ,PHAC(E)-PHAC(E)+180.;... 

FOR H-1 :9;„. 

N=ROUND(W(H)*2048/(20*2*PI)+1 )• 

JF REAL(FFTE(N)) < 0. ,PHAE(H)»PHAE(H)+180.;„. 

FOR G-1:9;... 

JFPHACfG) < 0. , PHAC(G)-PHAC(G)+360.;... 

FOr’k-1:9;... 

JF PHAE(K) < 0. . PHAE(K)«PHAE(K)+360.;„. 


PHAS = PHAC-PHAE; 
A=f.1 .0001,10 1;4.95 .250] 
PAGE 



WIND0W('211') 

//PLOT(A,'SCALE') 

PLOT(W,MDB, 'POINT-2*, *LOGLOG‘i 
YLABEL('MAGNITUDE’) } 

TITLE('FFT) 

WINDOW('212‘) 

PLOT('SCALE') 

PLOT(W, PHAS, 'POINT-2', ’LOG X') 
YLABEL('PHASE') } 

XLABEL('FREQUENCY RAD/SEC') 
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EllH(3 d0 2O47^' d6mi,iCa ' iOn USi " 9 m ° del 8> Z6r0 dela * with Was " 
E2=H(2:2046,2); 

E3*H(1 :2045,2); 

Cl=H(3:2047,3); 

C2=H(2:2046,3); 

C3-H(1:2045,3); 

C=C1;E=E1; 

V=H(4:2048,3); 

B«ONES(2045,1); 

A®[C1 ,C2,C3,E1 ,E2,E3,BJ; 

P*A\Y 

a 1 1 =p(1 , 1 ) ;a1 2-p(2, 1 );a1 3-p(3, 1 ) ; 
bl 1 =p(4, 1 ) ;b1 2-p(5, 1 );b1 3«p(6, 1 ) ; 

BIAS*(P(7,1)/(B1 1+B12+B13)) 

BB-BIAS*B; 

NU-[0B11 B12B13]; 

DE»[1 -All -A12-A13]; 

[A1 Bl Cl D1]-TF2SS(NU,DE); 

EE-E+BB; 

X»DSIM(A1 ,B1 ,C1 ,D1 ,EE’); 

XX-X*; 

R2-1 -(SUM((C-XX)‘*2))/SUM(C“2) 
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// •• 28d0b; LSE ide "" ,ic f us ' n 9 model 8 with a 3 time 
E1-H(3:2044,2); C ° nS,an( 

E2«H(2:2043,2); 

E3«H(1 :2042,2); 

C1»H(6:2047,3); 

C2=H(5:2046,3); 

C3=H(4:2045,3); 

C=Cl ;E=sEl ; 

Y=H(7:2048,3); 

MCl.C2.C3.E1.E2.E3J; 

p=a\v J ‘ 

Ri i *d! 1 ’ 1 ^ :A ^* P (2. 1 ) ;A1 3»P(3, 1 ); 

NU=[0 0 0 0 B1 1 B12 B131- 
DE=f1 -A1 1 -A12 -A13 0 0 0 1; 
v no 1 01 D1 J" TF2 SS(NU,DE)* 

X-DSIM(A1.B1,C1 I D1,E1* 

R2=1-((SUM((C-X')**2))/SUM((C**2))) 
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// 28d0b; LSE identif ''cation using model 8 

E1 “H( 3:2047.2); 

E2»H(2:2046,2); 

E3»H(1:2045,2); 

C1»H(3:2047,3); 

C2-H(2:2046,3); 

C3-H(1 :2045,3); 

C*C1;E-E1; 

V«H(4:2048,3); 

MC1.C2.C3.E1.E2.E3J; 

p*a\v J 

n 1 1 “a! 1 ' 1 A1 2 * p ( 2 > 1 ).’A1 3-P(3, 1 ) • 

NU=[0B11 B12B13J; ‘ 

DE»[1 -All -A12 -A13J; 

[A1 B1 Cl D1]-TF2SS(NU DEV 

X.DSlM(A1.Bl.Cl.Dl.En- 

R2-1-((SUM((C-X')**2))/sijM((C"2))) 



// z9d0, LSE identification using model 9 ** 

E1-H(4:2047,2); 

E2-H(3:2046,2); 

E3»H(2:2045,2); 

E4»H(1 :2044,2); 

CUH(4:2047,3); 

C2»H(3:2046,3); 

C3-H(2:2045,3); 

C*C1;E»E1; 

Y-H(5:2048,3); 

^■ C2 - C3 'E1.E2'E 3i E4J; 

n, 1 1 ’n! 1 ' 1);A12 - p (2.1);Al3-P(3,1); 

DE-[1 -All -A12-A13 0J; 

Vo 1 C1 Di J" TF 2 SS(NU,DE); 

X»DSIM(A1 ,B1 ,C1 ,D1 ,E’); 

R2«1-((SUM((C-X')**2))/SUM((C**2))) 



